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摘 要 : 天 文 软 件 开 发 与 应 用 中 迫切 需要 在 单机 环境 下 进行 高 性 能 的 科学 数据 处 理工 作 ,， 由 
于 机 器 配置 不 同 ， 采 用 传统 的 CUDA+GPU 技 术 存 在 明显 的 局 限 ， 不 利于 天 文 软 件 的 快速 移植 和 无 
颖 运行 。 针 对 明 安 图 频谱 射电 日 像 仪 (MingantU SpEctral Radioheliograph, MUSER) 数据 处 
理 中 的 网 格 化 (Gridding) 算法 ， 采 用 并 行 计 算 0penCL 技 术 进 行 多 线程 编程 实现 。 实 验 结果 表 
明 ， 基 于 0penCL 实 现 的 gridding 算 法 不 仅 能 够 在 多 种 GPU 上 运行 ， 而 且 能 够 在 纯 CPU 上 运行 。 当 
选择 在 GPU 上 执行 时 ， 算 法 的 执行 效率 与 基于 CUDA 实 现 的 gridding 算 法 执行 效率 大 致 相当 ,但 算 
法 不 局 限于 NVIDIA GPU, 解决 了 算法 对 CUDA+GPU 的 依赖 ; 同时 算法 也 能 在 纯 CPU 上 较 快 速 地 执行 ， 
适用 于 单机 模式 下 进行 天 文 软件 的 开发 和 测试 ， 也 便于 天 文 软件 的 应 用 与 推广 。 
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中 国 新 一 代 厘 米 - 分 米 波 射电 日 像 仪 一 明 安 图 频谱 射电 日 像 仪 MingantU SpEctral 
Radioheliograph, MUSER) 是 基于 利用 综合 孔径 成 像 原理 实现 在 频率 范围 0.4- 15. 0 GHz 内 对 太 
阳 进 行 高 分 辩 率 〈 时 间 、 空 间 ) 观测 并 成 像 的 射电 望远镜 一 ” 。 

MUSER 主要 由 天 线 ( 天 线 阵 ) 、 接 收 系统 、 数 据 处 理 系 统 三 部 分 组 成 。 天 线 阵 由 低频 阵 
(MUSER-I， 共 40 面 天 线 ) 和 高 频 阵 (MUSER-II， 共 60 面 天 线 ) 两 部 分 组 成 ， 天 线 阵 的 拓扑 结构 为 
ewe, E MUSER 观测 过 程 中 ， 低 频 阵 和 高 频 阵 分 别 每 3 ms 产生 一 个 数据 帧 ， 每 个 数据 帧 
分 别 用 100 KBytes 和 256 KBytes 存储 ， 数 据 流 量 分 别 为 32 MB/s 和 86 MB/s， 以 每 天 观测 10 
小 时 为 例 ， 每 日 原始 观测 数据 量 分 别 约 为 1.2 TB 和 3.2 TB”, 
如 此 海量 的 天 文 数据 ， 对 于 MUSER 数据 处 理工 作 而 言 是 一 个 巨大 的 挑战 ， 但 挑战 与 机 遇 并 
存 。 在 前 期 研究 工作 中 ， 为 了 满足 数据 处 理 需 求 ， 研 究 人 员 一 方面 采用 分 布 式 计算 技术 ， 设 计 
并 实现 了 一 套 分 布 式 数据 处 理 执 行 框架 OpenCluster”; 另 一 方面 采用 并 行 计算 技术 并 借助 高 性 
能 计算 设备 ， 专 注 于 底层 算法 研究 ， 实 现 了 一 套 高 效 的 射电 于 涉 阵 数据 处 理 管线 ”。 特 别 地 ， 
采用 MPI (Message Passing Interface， 消息 传递 接口 ) CUDA (Compute Unified Device 
Architecture， 统 一 计算 设备 架构 ) 以 及 OpenCL (Open Computing Language, a 
技术 已 经 在 天 文中 有 相关 的 应 用 ， 并 已 经 取得 一 定 的 成 果 ”，…"。 证 明 这 些 技术 稳定 性 与 性 能 

以 满足 天 文 应 用 软件 开发 的 要 求 。 

为 了 便于 MUSER 数据 处 理 系统 的 开发 、 部 署 、 应 用 以 及 推广 ， 同 时 满足 单机 环境 下 仍旧 能 
够 进行 高 性 能 的 成 像 数据 处 理工 作 。 在 前 期 对 OpenCL 技术 研究 基础 上 ， 进 一 步 采用 OpenCL 技 
术 对 射电 干涉 阵 成 像 中 的 网 格 化 (Gridding) 算法 进行 并 行 优 化 ,在 保证 gridding 算法 执行 效 
率 的 同时 ， 提 升 算法 对 各 种 硬件 平台 的 适应 性 ， 便 于 后 期 MUSER 数据 处 理 软件 的 应 用 与 推广 。 
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1 基于 0penCL 的 gridding 算 法 实现 


1.1 Gridding 算法 
射 
(Gridding), 


E F H E RARA TRAS 


等 操作 ， 最 后 产生 射 


为 了 借助 离散 侍 
算 效率 ， 进 行 快速 成 像 ， 射 
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孔径 成 像 原 理 对 采样 


图 1 左 图 所 示 )。 实 际 


包干 涉 阵 采集 


的 可 见 度数 据 进 行 网 格 化 


快速 傅 里 时 变换 (Fast Fourier Transformation, FFT)、 去 卷 积 (Deconvolution) 
电源 的 观测 图 像 。 

里 叶 变 换 (Discrete Fourier Transform，DFT) 的 快速 方法 一 FFT 一 的 计 
到 的 可 见 度数 据 必 须 落 在 一 个 均匀 划分 
中 ， 由 于 UV 采样 点 分 布 不 均匀 ， 干 涉 仪 得 到 的 是 非 均 匀 采 样 可 见 度 数据 


的 网 格 上 《如 


(Non-uniform Sampled Visibility Data)〔( 如 图 1 右 图 所 示 )， 将 非 均匀 采样 可 见 度数 据 插值 
到 均匀 划分 的 网 格 点 上 的 过 程 称 为 网 格 化 Gridding) 099%, 


Fig. 1 Regular grid (left) and the schema 


将 非 均匀 采样 可 见 度数 据 转 换 成 均匀 采样 可 见 度 数据 ,有 


插值 ， 双 线性 插值 等 


积 核 法 ， 出 于 实际 考虑 ， 


主要 有 两 类 方法 : 内 插 法 和 卷 积 核 法 。 


一 些 内 插值 方法 被 开发 


; 文献 [15] 表 明 最 优 的 网 格 化 方法 是 卷 积 


有 限 的 卷 积 函数 取代 了 无 限 扩展 的 
(Gridding Convolution Function, GCF) 有 截断 sinc 


ic diagram of non-uniform visibi 


对 未 知 像 


CA) 


= 1 


en > 


图 1 均匀 划分 的 网 格 〈 左 ) 与 非 均 匀 分 布 的 可 见 度数 据 示 意图 


ity data (right) 


数 ， 文 献 [16] 讨 论 了 这 些 网 格 卷 积 函 数 在 不 同 尺寸 大 小 下 的 1 
选取 了 大 小 为 6 的 拟 合 的 球面 波 函 数 作为 网 格 卷 积 函 数 。 


/ Non-uniform 
/ sampled 
Á visibility data 


1.2 Gridding 算法 并 行 化 
程 规范 ，0penCL FE 


根据 OpenCL 9 


GCF 


Fig. 2 Conv 
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素 点 的 数据 进行 重 采 样 ， 
于 重建 非 均 匀 采 样 ， 如 最 邻近 
一 个 无 限 扩展 的 sinc 函数 ， 即 卷 
sinc PRX. FA AY PYAR EAR ER EL 
函数 ， 截 断 指数 函数 ， 
“EAE. ZE MUSER 数据 处 理 


拟 合 的 球面 波 函 
中 ， 我 们 


Regular uv- 
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grid data J 
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Dirty / 


7 image 


12 卷 积 网 格 化 并 成 像 


utional gridding and ima 


这 主 要 由 两 部 分 组 成 : 


ging 
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部 分 是 


AN 


用 


OpenCL 语言 (类 CC 


语言 ) 编写 的 内 核 函数 ， 执 行 在 OpenCL 设备 (多 核 CPU、GPU、Cell、DSP、FPGA 等 ) E; A— 


部 分 是 通过 调用 OpenCL 的 应 用 编程 接口 (Application Programing Interface, API) # 
了 在 主机 (CPU) 上 。 


核 函 数 的 主机 程序 ， 


Há 


本 文 将 gridding 过 程 中 的 关键 步骤 ,由 OpenCL 语言 编写 成 内 核 函 数 ， 


管理 内 


以 便 在 OpenCL 设备 
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上 并 行 执行 ， 再 通过 在 主机 上 调用 这 些 内 核 函 数 ， 从 而 实现 gridding 算法 的 并 行 化 。 主 要 改写 
成 的 内 核 函 数 有 : gridding kernel , correct grid kernel , normalize kernel , 
point symmetric_kernel, shift kernel 等 。 此 外 ， 成 像 过 程 的 傅 里 时 变换 是 通过 调用 基于 


OpenCL 实现 的 快速 全 里 叶 变 换 FFT， 包 含 于 python 的 pyfft 包 中 。 
GCFAGE 


非 均匀 采样 ; 
可 见 度数 据 均匀 UV 网 格 数据 


读 取 UVFITS 文 件 


去 着 各 脏 图 及 脏 束 


图 3 Gridding 并 行 执行 流程 图 
Fig.3 The flow diagram of parallel gridding 


本 文 处 理 的 图 和 像 数 据 是 以 二 维 数组 的 形式 表示 ， 在 采用 0penCL 语言 编写 内 核 函 数 过 程 中， 

由 OpenCL 提供 的 二 维 数组 索引 get global id(0) 和 get global id(1)， 唯 一 标识 二 维 数 组 中 

的 元 素 ， 为 方便 计算 ， 需 借助 传递 给 内 核 函 数 的 二 维 数组 的 参数 〈 宽 度 或 高 度 )， 将 二 维 数组 按 

行 或 按 列 转换 成 一 维 数组 ,这样 ， 二 维 索引 get lebal id 和 set -8lebal id 人 转换 成 一 维 

索引 逢 id， 当 内 核 函 数 在 并 行 设备 中 执行 时 ， 每 个 线程 处 理由 索引 和 值 -id 唯一 标识 的 二 维 数组 
元 素 。 


id=ix*W+iy 
=3*442 


图 4 OpenCL 二 维 索引 转换 一 维 索引 示意 图 


Fig. 4 The schematic diagram of OpenCL two-dimensional index to one-dimensional index 
在 gridding 执行 过 程 中 , E FE AAA AER AA ARG, 
为 减少 网 格 化 的 计算 量 ， 先 对 一 半 的 可 见 度数 据 进行 插值 网 格 化 ， 再 通过 中 心 点 对 称 方式 计算 
得 到 男 一 半数 据 。 将 中 心 点 对 称 操作 编写 成 内 核 函 数 point_symmetric kernel， 示 意图 如 下 。 


iy < W/2 


ix=get global id(0) 


0+0i 0+0i 3+3i 4+4i 0+0i 0+0i 3+3i 4+4i 
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图 5 中 心 点 对 称 操作 示意 图 
Fig.5 The schematic diagram of point symmetric operation 


point symmetric kernel 伪 代 码 如 下 : 


OpenCL kernel function: point symmetric kerne 


Input parameter: a two dimensional array im, a flag hfac, the two dimensional array height H and width W 
//Obtain two dimensional array index (ix, iy) 
ixeget global id(0) 
iyeget global_id(1) 
// Stay within the bounds 
If ix > 0 and ix < H and iy > 0 and iy < W/2 then 
nix«H-ix 
niy<W-iy 
im[ix*W + iy]. x—im[nix*W + niy]. x 
im[ix*W + iy]. y—hfacxim[nix*W + niy]. y 
Output: a new two dimensional array im 


在 傅 里 叶 变 换 过 程 中 , 为 使 变换 后 的 图 像 是 一 个 完整 的 周期 ， 需 要 进行 平移 (shift ) 操作 。 
将 shift 操作 编写 成 内 核 函 数 shift kernel1， 示 意图 如 下 。 


A AR 


1+i < 2+2i 3+3i 4+4i | 11+11i 12+12i 9+9i 10+10i 


5+5i 6+6i 7+7i 8+8i i 15+15i 16+16i 13+13i 14+14i 
shift_kernel 
9+9i 10+10i 11+11i 12+12i 3+3i 4+4i 1+1i 2+2i 
13+13i 14+14i 15+15i 16+16i 7+Ti 8+8i 5+5i 6+6i 


图 6 傅 里 叶 变 换 平移 操作 示意 图 


Fig.6 The schematic diagram of Fourier transform shift operation 
在 gridding 执行 过 程 中 ， 由 于 可 见 度 数据 与 一 个 网 格 卷 积 函 数 GCF 进行 了 卷 积 ， 为 消除 
GCF 的 影响 ， 需 进行 网 格 校正 操作 。 将 网 格 校正 操作 编写 成 内 核 函 数 correct_grid kernel。 些 
外 一 进行 网 格 校正 后 ， 需 要 对 脏 图 和 脏 束 进 一 步 的 进行 归 一 化 或 标准 化 一 将 脏 束 进行 思 一 化 ， 
即将 脏 图 和 脏 束 中 所 有 元 素 除 以 脏 束 中 的 最 大 值 一 对 脏 图 进行 标准 化 一 即 以 及 将 脏 图 中 所 有 所 
以 元 素 除 以 脏 束 中 的 最 夫 值 ， 从 而 保证 脏 束 的 峰值 为 1 以 及 脏 图 与 脏 束 具有 相同 量 纲 标准 化 后 
的 脏 束 和 脏 图 具有 相同 的 量 纲 ， 便 于 后 续 进 行 洁 化 。 


2 4 5 2 1snum 4*num S*num enum 
| 
和 一 | 
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normalize kernel 1*num 3*num 4*num 1snum 
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rT | i T T 


1.3 并 行 加 权 

由 于 UV 采样 点 太 少 ， 对 可 见 度 函数 采样 值 进行 傅 里 叶 首 变换 得 到 的 图 像 ( 称 汶 脏 图 ) 包含 
一 些 虚 假 信息 , 通过 对 UV 采样 点 赋予 不 同 的 权 值 来 改善 图 像 质量 的 操作 称 之 为 加 权 ” 。 常 用 的 
加 权 方 式 有 自然 加 权 , 统一 加 权 以 及 Taper HAL, 本 文 将 这 三 种 加 权 方 式 改 写成 OpenCL 内 核 函 
数 weight_natural kernel, weight uniform kernel Fl weight taper kernel. 


三 种 加 权 方 法 伪 代 码 如 下 : 
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OpenCL kernel function: weight natural kernel、 weight uniform kernel、 weight taper kernel 


Input parameter: a two dimensional array nim, the number of samples falling into cell cnt, th 
height H and width W 
//Obtain two dimensional array index (ix, iy) 
ixeget_global_id(0) 
iyeget global id(1) 
id—iy + ix*W 
// Stay within the bounds 
If ix > 0 and ix < H and iy > 0 and iy < W then 
Af entlid] 15 0 then 0 
//natural weighting 
wet<l. 
//uniform weighting 
wgte—1./ entlid] 
//taper weighting 
H wet—exp (- ((ix-H/2) * (ix-H/2) + (iy-W/2) * (iy-W/2)) / (2x0. 5*0. 5*WxW) ) 
aid 
nimlid]. y—nimLid]. y*wgt 
Output: weighted nim 


e two dimensional array 


自然 加 权 赋 予 相 应 像素 点 的 权重 为 1, RIMA RR ANE 


为 区 域内 采样 点 频数 


的 倒数 ，Taper 加 权 赋 予 相 应 像素 点 的 权重 为 对 应 高 斯 函数 值 。 通 过 改变 像素 点 的 权重 大 小 ， 对 


可 见 度 数据 进行 加 权 操 作 ， 从 而 改善 图 像 质量 


o 


2 实验 结果 


2.1 实验 环境 


OpenCL 定义 不 同 的 公司 或 厂商 为 不 同 的 平台 Platform) ， 每 个 公 


司 或 厂商 提供 不 同 的 


OpenCL 设备 (Device) 。 本 文选 用 的 OpenCL 设备 有 NVIDIA 的 GPU (GeForce GTX TITAN X, Tesla 
k20m) 和 Intel 的 CPU (Intel Xeon E5-2620 V2) ， 平 台 及 设备 参数 见 表 2. 


表 2 平台 和 设备 参数 信息 


Tab.2 The information of platform and device 


Parameter Intel CPU GTX TITAN X GPU Tesla k20m GPU 
Platform vendor Intel (R) Corporation NVIDIA Corporation NVIDIA Corporation 
latform version OpenCL 1.2 LINUX OpenCL 1.2 CUDA 7.5.18 OpenCL 1.2 CUDA 7.5.18 
evice name Intel (R) Xeon(R) CPU E5-2620 GeForce GTX TITAN X Tesla k20m 

V2 @ 2. 10GHz 
evice type CPU GPU GPU 
evice max clock speed 2100 MHz 1215 MHz 705 MHz 
evice compute units 24 24 13 
evice cores 6 2496 2496 
evice max work group size 8192 1024 1024 


2.2 实验 结果 
实验 数据 来 源 于 MUSER 2015 E11 A 1 H 12H) 8 4) 494) 354 毫秒 时 
原始 观测 数据 经 过 一 系列 预 处 理 〈 坏 天 线 标记 ， 星 表 碍 询 等 ) 合成 为 标准 


刻 对 太阳 的 观测 数据 ， 
天 文 数据 格式 (UVFITS 


文件 ) ， 数 据 文 件 被 命名 为 20151101-120849 354161240. uvfits。 通 过 MUSER 数据 处 理 系统 执 
行 网 格 化 和 快速 傅 里 叶 变 换 等 操作 ， 生 成 相应 时 刻 太阳 亮度 分 布 图 像 ( 未 经 过 去 卷 积 操作 的 图 
像 称 为 脏 图 ) 以 及 对 应 的 点 扩散 函数 (Point Spread Function, PSF), PSF 也 称 为 脏 束 。 


图 87 MUSER 成 的 脏 图 〈 左 和 点 扩散 函数 PSF〈 右 )， 观 测 时 间 是 2015 4 


TIME: 2015-11-01T12:08:49.354161240 POL: R @1.7125GHz 
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E11 月 1 日 12:08:49:354， 频 率 是 1.7125 GHz， 极 化 方式 


Fig.87 The dirty image (left) and PSF (right) of MUSER at 12:08:49:354 on November 11th, 2015 (Frequency: 1.7125 GHz, 


Polarization: right) 


在 同一 台 服 务 器 上 ， 分 别 在 CPU (Intel Xeon E5-2620 V2), GPU (GeForce GTX TITAN X, 


Tesla k20m) 并 行 设 备 上 执行 基于 OpenCL 实现 的 gridding 算法 和 在 GPU (GeForce GTX TITAN X, 


Tesla k20m) 上 执行 基于 CUDA 实现 的 gridding 算法 ， 并 进行 快速 傅 里 时 变换 成 图 ， 分 别 生成 大 


小 为 1024*1024 pixels 和 2048*2048 pixels 的 及 


行 时 间 参 照 下 表 。 


# 3 Gridding+IFFT 过 程 执行 时 间 表 
Tab.3 the Execution time of Gridding and FFT 


图 。 网 格 化 并 成 图 (Gridding+IFFT) 过 程 执 


Device 设备 Image size 图像 大 小 Bxeetutien time-Cs》 执 行 时 间 (s) 
CPU OpenCL 1024*1024 0. 904 
2048*2048 4.407 
GPU OpenCL (TITAN X) 1024*1024 0. 350 
2048*2048 0.713 
GPU OpenCL (k20m) 10241024 0. 372 
20482048 0. 782 
GPU CUDA (TITAN X) 1024*1024 0. 323 
20482048 0.610 
GPU CUDA (k20m) 1024*1024 0. 344 
2048*2048 0. 760 
从 表 中 可 以 看 出 , 执行 基于 OpenCL 和 基于 CUDA 实现 的 gridding 算法 并 成 图 , 在 GPU 环境 
下 所 消耗 的 时 间 基 本 相当 ， 同 时 ， 基 于 OpenCL 实现 的 gridding 算法 还 能 在 纯 CPU 环境 下 较 快 
速 地 被 执行 。 
2.3 讨论 
本 文采 用 并 行 计算 OpenCL BER, 基于 Python 语言 导入 PyOpenCL 与 PyFFT 扩展 包 , 对 MUSER 
成 像 过 程 中 的 gridding 算法 以 及 快速 传 里 叶 变换 FFT 成 图 过 程 进行 了 实现 。 基 于 OpenCL 实现 


的 gridding 算法 
(1) 能 够 在 GPU 环境 中 执行 ， 并 且 不 局 限 与 NVIDIA GPU (实验 条 件 限 秆 


只 有 如 下 优点 : 


= 


J NVIDIA 


GPU 和 Intel CPU), 执行 效率 与 基于 CUDA 实现 的 gridding 算法 执行 效率 大 致 相当 , 保 


成 图 的 性 能 ; 


> REH 


UE 了 MUSER 


(2) 能 够 在 多 核 CPU 环境 下 执行 ， 适 用 于 在 单机 环境 下 进行 开发 与 测试 ， 便 于 MUSER 软件 


的 推广 与 应 


Jo 


综 上 所 述 ， 采 用 OpenCL 实现 的 并 行 gridding 算法 ， 在 保证 算法 执行 效率 的 同时 有 效 地 提 
升 了 算法 对 硬件 平台 的 适应 性 。 此 外 , 本 工作 进一步 验证 了 0penCL 在 天 文 软件 开发 中 的 可 用 性 ， 
从 CUDA 到 OpenCL， 异 构 系 统 从 NVIDIA GPU+CPU 的 异 构 模 式 转变 成 并 行 设备 +CPU 的 异 构 模 式 。 


这 种 异 构 模 式 的 转变 将 扩展 并 行 计算 在 天 文 高 性 能 软件 开发 中 的 应 用 


3 下 一 步 工 作 


本 文 研究 的 射电 干涉 阵 成 像 过 程 局 限于 小 视 场 ， 在 大 视 场 成 像 中 w 项 引起 的 视 场 扭曲 是 不 
可 忽略 的 ， 传 统 的 二 维 FFT 成 像 将 不 再 适用 ， 因 此 ， 下 一 步 我 们 将 在 此 研究 工作 的 基础 上 ， 采 
用 OpenCL 技术 对 大 视 场 成 像 中 的 关键 算法 ， 例 如 w-projection，faceting 以 及 大 视 场 混合 
去 ， 进 行 并 行 优化 。 
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Abstract: It's urgent to carry out high-performance scientific data processing with a single 
machine in the development and application of astronomical software. However, due to the different 
configurations of machines, the traditional CUDA + GPU technology has obvious limitations in 
portability and seamlessness. According to gridding algorithm in MingantU SpEctral Radioheliograph 
(MUSER) data processing, the OpenCL technology is used in parallel to implement multi-thread 
programming. The experimental results show that the gridding algorithm based on OpenCL can not 
only can run on varieus-GPUs, but also merely on CPUs. While choosing execution on GPU, the 
execution efficiency of gridding algorithm is approximately equal with 1t based on CUDA. At the 
same time, the algorithm is not limited to the NVIDIA’s GPU, which has solved the problem of 
environmental dependence of CUDA+GPU. And the algorithm also has an acceptable execution 
efficiency high+mplementation with the merely CPU, which is suitable for development and testing 
astronomy software with a single machine —but-alse and will facilitate the application and promotion 
of astronomical software. 

Key words: Gridding; Parallel computing; OpenCL; MUSER 
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